// Test for returning 4momentum clonearray
#include "flow_do.h"
#include "FlowToyMCBase.h"
#include "TH1D.h"
#include <vector>
#include <iostream>
#include "TROOT.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TF1.h"
#include "TFile.h"
#include "TChain.h"
#include "TClonesArray.h"
#include "TLorentzVector.h"
#include <TStopwatch.h>
#include "TRandom3.h"

using namespace std; // 




int main(int argc, char **argv)
{

		TROOT root("flow","mc");
		if ( argc<6 ) {
				cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
				cout<<"usage :  "<<argv[0]<<" <outputFile> <N_target> <N_EP> <Rndm_seed> <Nevt>"<<endl;
				cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
				cout << endl << endl;
				exit(1);
		}

		// CONSTANT
		char *outFile = argv[1];
		int N_ch[ND];
		N_ch[kTargetDetector] = atoi(argv[2]);  
		int N_EPdetctor = atoi(argv[3]); 
		int Rndm_seed = atoi(argv[4]);
		int numberEvents= atoi(argv[5]);

		for(int i= 1; i<ND; i++){
 			N_ch[i] = N_EPdetctor;
		}

		TFile *fout = new TFile(outFile,"recreate");

		cout<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
		cout << "Outfile : "<< outFile << endl;
		cout << "Inputs : "<<endl;
		cout <<"\t nPaarticle for Target   = "<< N_ch[0] << endl;
		cout <<"\t nPaarticle for EP Source= "<< N_ch[1] << endl;
		cout <<"\t nEvt = "<< numberEvents  << endl;
		cout<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;

		// ---------- prepareing histograms --------------
		TH1D *hphidist[ND];
		TH1D *hpsidist[NH];
		TH1D *hpsi_caldist[ND][NH];
		TH1D *hdpsi[ND][NH];
		TH1D *hResol[ND][NH];
		TH1D *hVn[ND][NH];
		TH1D *hVn1evt[ND][NH];
		TH1D *hVnReal[ND][NH];
		TH1D *hphidist_measured[ND];
		TH1D *hE[ND];
		TH1D *hPt[ND];
		TH1D *hSub2Resol[ND][NH];
		TH1D *hSub3Resol[ND][NH][3];
		TRandom3 *MyRandom;

		double EP2[ND]={};
		double EP3[ND]={};
		double* Sub3Resol[ND][NH];

		for(int ih=0; ih<4; ih++){
			hpsidist[ih] = new TH1D(Form("hpsidist%02d", ih), "" , 1024, 0, 2*pi);
		}
		for (int id=0; id<ND; id++) {	
				hphidist[id] = new TH1D(Form("hphidist%02d",id),"phi distribution over all event", 1024, 0, 2*pi);
				hphidist_measured[id] = new TH1D(Form("hphidist_measured%02d",id),"phi_measured distribution over all event", 1024, 0, 2*pi);
				hE[id] = new TH1D(Form("hE%02d",id), "",1024,0,10);			
				hPt[id] = new TH1D(Form("hPt%02d",id),"",1024,0,20);	
			for(int ih=0; ih<4; ih++){

				hpsi_caldist[id][ih] = new TH1D(Form("hpsi_caldist%02d_%02d",id, ih), "" , 1024, -2*pi, 2*pi);
				hdpsi[id][ih] = new TH1D(Form("hdpsi%02d_%02d", id, ih), "" , 1024, -2*pi, 2*pi);
				hResol[id][ih] = new TH1D(Form("hResol%02d_%02d",id, ih), "" , 1024, -1, 1);
				hVn[id][ih] = new TH1D(Form("hVn%02d_%02d",id, ih), "" , 1024, -3, 3);	
				hVn1evt[id][ih] = new TH1D(Form("hVn1evt%02d_%02d",id, ih), "" , 1024, -3, 3);	
				hVnReal[id][ih] = new TH1D(Form("hVnReal%02d_%02d",id, ih), "" , 1024, -3, 3);
				hSub2Resol[id][ih] = new TH1D(Form("hSub2Resol%02d_%02d",id,ih),"",1024,-1,1);
				for(int isub=0;isub<3;isub++){
					hSub3Resol[id][ih][isub] = new TH1D(Form("hSub3Resol%02d_%02d_%02d",id,ih,isub),"",1024,-1,1);
					}
			}
		}
		MyRandom = new TRandom3();

		TStopwatch timer;
		timer.Start();
		FlowToyMCBase *mc = new FlowToyMCBase();
		mc->SetEPSeed(Rndm_seed);
		//*****************************************************************************************//
		// ---------------------  event loop start --------------- //
		//*****************************************************************************************//
		for(int ievt =0; ievt<numberEvents; ievt++){
				if( ievt%100==0) {
						cout << " Event: " << ievt << " " << ievt*100.0/numberEvents <<"% "<< endl;
				}
				mc->SetEP();
				for (int id =1; id<ND; id++) {
					mc->Init();
					mc->SetNPar(N_ch[id]);
					mc->SetDetResol(det_resol[id]);						
					mc->Generate4Vector();
					for (int ip=0; ip<N_ch[id];ip++) {
						hphidist_measured[id]->Fill(mc->GetPhi(ip));
					}
					EP2[id] = mc->Cal_EP(2);		
					hpsi_caldist[id][2]->Fill( EP2[id] );
					EP3[id] = mc->Cal_EP(3);
					hpsi_caldist[id][3]->Fill( EP3[id] );
					hSub2Resol[id][2]->Fill( mc->Cal_Sub2Resol(2) );
					hSub2Resol[id][3]->Fill( mc->Cal_Sub2Resol(3) );
					for(int isub=0;isub<3;isub++){
						int sub1 = MyRandom->Uniform(0,N_ch[id]);
						int sub2 = MyRandom->Uniform(0,N_ch[id]);
						while (sub1 + sub2 >= N_ch[id]){
							sub1= MyRandom->Uniform(0,N_ch[id]);
							sub2= MyRandom->Uniform(0,N_ch[id]);
							}
						
						Sub3Resol[id][2]=mc->Cal_Sub3Resol(2, sub1, sub2);///////////////	
						hSub3Resol[id][2][isub]->Fill(Sub3Resol[id][2][isub]);
						Sub3Resol[id][3]=mc->Cal_Sub3Resol(3, sub1, sub2);///////////////	
						hSub3Resol[id][3][isub]->Fill(Sub3Resol[id][3][isub]);
						}
					hResol[id][2]->Fill( mc->Cal_Resol(2) );	
					hResol[id][3]->Fill( mc->Cal_Resol(3) );
					for(int ip=0;ip<mc->GetNPar();ip++){
					hE[id]->Fill(mc->GetE(ip));
					hPt[id]->Fill(mc->GetPt(ip));
					}	
				}//detector loop
					mc->Init();
					mc->SetNPar(N_ch[0]);
					mc->SetDetResol(0); // Ideal
					mc->GeneratePhi();
				for (int id=1; id<ND; id++){ // Target detector
					mc->SetEP_Measured(2, EP2[id]);
					mc->SetEP_Measured(3, EP3[id]);
					hVn[id][2]->Fill( mc->Cal_Vn(2) );	
					hVn[id][3]->Fill( mc->Cal_Vn(3) );
					hdpsi[id][2]->Fill( mc->GetEP_Measured(2) - mc->GetEP_Input(2) );
					hdpsi[id][3]->Fill( mc->GetEP_Measured(3) - mc->GetEP_Input(3) );
					hVnReal[id][2]->Fill( mc->Cal_VnReal(2) );	
					hVnReal[id][3]->Fill( mc->Cal_VnReal(3) );
				}
				// Generation done. ( 0: target(TPC), 1: source(For EP measure) )	


		} // end of evnetp loop

		// ***********************************************
		// Save Histo ***********************************
		// **********************************************
		for(int id=1; id<ND; id++){
				hphidist_measured[id]->Write();	
				hE[id]->Write();
				hPt[id]->Write();
				for( int ih=2; ih < NH ; ih++ ){
						hpsi_caldist[id][ih]->Write();
						hdpsi[id][ih]->Write();
						hResol[id][ih]->Write();
						hVn[id][ih]->Write();
						hVnReal[id][ih]->Write();	
						hSub2Resol[id][ih]->Write();
						for(int isub=0;isub<3;isub++){
							hSub3Resol[id][ih][isub]->Write();
							}
				}			
		} 
		cout <<"Total number of event processed = "<< numberEvents << endl;
		cout <<"Writing histograms into a root file"<<endl;
		fout->cd();
		fout->Close();
		timer.Print();
		cout <<"Sucessfully finished."<<endl;
}
